home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Source Code / Libraries / DCLAP 4j / SeqPups / apps / clustalw.src / sequence.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-12-17  |  14.3 KB  |  597 lines  |  [TEXT/R*ch]
  1. /********* Sequence input routines for CLUSTAL W *******************/
  2. /* DES was here.  FEB. 1994 */
  3. /* Now reads PILEUP/MSF and CLUSTAL alignment files */
  4.  
  5. #include <stdio.h>
  6. #include <string.h>
  7. #include <ctype.h>
  8. #include <stdlib.h>
  9. #include "clustalw.h"    
  10.  
  11. #define MIN(a,b) ((a)<(b)?(a):(b))
  12.  
  13. /*
  14. *    Prototypes
  15. */
  16.  
  17. extern Boolean linetype(char *,char *);
  18. extern Boolean blankline(char *);
  19. extern void warning(char *,...);
  20. extern void error(char *,...);
  21. extern char *    rtrim(char *);
  22. extern char *    blank_to_(char *);
  23. extern void     getstr(char *,char *);
  24.  
  25. void fill_chartab(void);
  26. int readseqs(int);
  27.  
  28. static void         get_seq(char *,char *,int *,char *);
  29. static void get_clustal_seq(char *,char *,int *,char *,int);
  30. static void     get_msf_seq(char *,char *,int *,char *,int);
  31. static void check_infile(int *);
  32. static void p_encode(char *, char *, int);
  33. static void n_encode(char *, char *, int);
  34. static int res_index(char *,char);
  35. static Boolean check_dnaflag(char *, int);
  36. static int count_clustal_seqs(void);
  37. static int count_msf_seqs(void);
  38.  
  39. /*
  40.  *    Global variables
  41.  */
  42.  
  43. extern FILE *fin;
  44. extern Boolean usemenu, dnaflag, explicit_dnaflag;
  45. extern char seqname[];
  46. extern int nseqs;
  47. extern int *seqlen_array;
  48. extern int *output_index;
  49. extern char **names,**titles;
  50. extern char **seq_array;
  51. extern profile1_empty, profile2_empty;
  52. extern int max_aa;
  53. extern int gap_pos2;
  54.  
  55. char *amino_acid_codes   =    "ABCDEFGHIKLMNPQRSTUVWXYZ-";  /* DES */
  56. char *nucleic_acid_order =       "ACGTUN";
  57. static int seqFormat;
  58. static char chartab[128];
  59. static char *formatNames[] = {"unknown","EMBL/Swiss-Prot","PIR",
  60.                   "Pearson","GDE","Clustal","Pileup/MSF","USER"};
  61.  
  62. void fill_chartab(void)    /* Create translation and check table */
  63. {
  64.     register int i;
  65.     register char c;
  66.     
  67.     for(i=0;i<128;chartab[i++]=0);
  68.     for(i=0;c=amino_acid_codes[i];i++)
  69.         chartab[c]=chartab[tolower(c)]=c;
  70. }
  71.  
  72. static void get_msf_seq(char *sname,char *seq,int *len,char *tit,int seqno)
  73. /* read the seqno_th. sequence from a PILEUP multiple alignment file */
  74. {
  75.     static char line[MAXLINE+1];
  76.     int i,j,k;
  77.     unsigned char c;
  78.  
  79.     fseek(fin,0,0);         /* start at the beginning */
  80.  
  81.     *len=0;                /* initialise length to zero */
  82.         for(i=0;;i++) {
  83.         if(fgets(line,MAXLINE+1,fin)==NULL) return; /* read the title*/
  84.         if(linetype(line,"//") ) break;            /* lines...ignore*/
  85.     }
  86.  
  87.     while (fgets(line,MAXLINE+1,fin) != NULL) {
  88.         if(!blankline(line)) {
  89.  
  90.             for(i=1;i<seqno;i++) fgets(line,MAXLINE+1,fin);
  91.                         for(j=0;j<=strlen(line);j++) if(line[j] != ' ') break;
  92.             for(k=j;k<=strlen(line);k++) if(line[k] == ' ') break;
  93.             strncpy(sname,line+j,MIN(MAXNAMES,k-j)); 
  94.             sname[MIN(MAXNAMES,k-j)]=EOS;
  95.             rtrim(sname);
  96.                            blank_to_(sname);
  97.  
  98.             for(i=k;*len < MAXLEN;i++) {
  99.                 c=line[i];
  100.                 if(c == '.') c = '-';
  101.                 if(c == '*') c = 'X';
  102.                 if(c == '\n' || c == EOS) break; /* EOL */
  103.                 if( (c=chartab[c])) seq[++(*len)]=c;
  104.             }
  105.             if(*len == MAXLEN) return;
  106.  
  107.             for(i=0;;i++) {
  108.                 if(fgets(line,MAXLINE+1,fin)==NULL) return;
  109.                 if(blankline(line)) break;
  110.             }
  111.         }
  112.     }
  113. }
  114.  
  115.  
  116. static void get_clustal_seq(char *sname,char *seq,int *len,char *tit,int seqno)
  117. /* read the seqno_th. sequence from a clustal multiple alignment file */
  118. {
  119.     static char line[MAXLINE+1];
  120.     int i,j,k;
  121.     unsigned char c;
  122.  
  123.     fseek(fin,0,0);         /* start at the beginning */
  124.  
  125.     *len=0;                /* initialise length to zero */
  126.     fgets(line,MAXLINE+1,fin);    /* read the title line...ignore it */
  127.  
  128.     while (fgets(line,MAXLINE+1,fin) != NULL) {
  129.         if(!blankline(line)) {
  130.  
  131.             for(i=1;i<seqno;i++) fgets(line,MAXLINE+1,fin);
  132.                         for(j=0;j<=strlen(line);j++) if(line[j] != ' ') break;
  133.             strncpy(sname,line+j,MAXNAMES-j); /* remember entryname */
  134.             sname[MAXNAMES]=EOS;
  135.             rtrim(sname);
  136.                            blank_to_(sname);
  137.  
  138.             for(i=14;*len < MAXLEN;i++) {
  139.                 c=line[i];
  140.                 if(c == '\n' || c == EOS) break; /* EOL */
  141.                 if( (c=chartab[c])) seq[++(*len)]=c;
  142.             }
  143.             if(*len == MAXLEN) return;
  144.  
  145.             for(i=0;;i++) {
  146.                 if(fgets(line,MAXLINE+1,fin)==NULL) return;
  147.                 if(blankline(line)) break;
  148.             }
  149.         }
  150.     }
  151. }
  152.  
  153.  
  154. static void get_seq(char *sname,char *seq,int *len,char *tit)
  155. {
  156.     static char line[MAXLINE+1];
  157.     int i, offset;
  158.         unsigned char c;
  159.  
  160.     switch(seqFormat) {
  161.  
  162. /************************************/
  163.         case EMBLSWISS:
  164.             while( !linetype(line,"ID") )
  165.                 fgets(line,MAXLINE+1,fin);
  166.             
  167.                         for(i=5;i<=strlen(line);i++)  /* DES */
  168.                 if(line[i] != ' ') break;
  169.             strncpy(sname,line+i,MAXNAMES); /* remember entryname */
  170.             sname[MAXNAMES]=EOS;
  171.             rtrim(sname);
  172.                         blank_to_(sname);
  173.  
  174.             while( !linetype(line,"SQ") )
  175.                 fgets(line,MAXLINE+1,fin);
  176.             
  177.             *len=0;
  178.             while(fgets(line,MAXLINE+1,fin)) {
  179.                 for(i=0;*len < MAXLEN;i++) {
  180.                     c=line[i];
  181.                 if(c == '\n' || c == EOS || c == '/')
  182.                     break;            /* EOL */
  183.                 if( (c=chartab[c]))
  184.                     seq[++(*len)]=c;
  185.                 }
  186.             if(*len == MAXLEN || c == '/') break;
  187.             }
  188.         break;
  189.         
  190. /************************************/
  191.         case PIR:
  192.             while(*line != '>')
  193.                 fgets(line,MAXLINE+1,fin);            
  194.                         for(i=4;i<=strlen(line);i++)  /* DES */
  195.                 if(line[i] != ' ') break;
  196.             strncpy(sname,line+i,MAXNAMES); /* remember entryname */
  197.             sname[MAXNAMES]=EOS;
  198.             rtrim(sname);
  199.                         blank_to_(sname);
  200.  
  201.             fgets(line,MAXLINE+1,fin);
  202.             strncpy(tit,line,MAXTITLES);
  203.             tit[MAXTITLES]=EOS;
  204.             i=strlen(tit);
  205.             if(tit[i-1]=='\n') tit[i-1]=EOS;
  206.             
  207.             *len=0;
  208.             while(fgets(line,MAXLINE+1,fin)) {
  209.                 for(i=0;*len < MAXLEN;i++) {
  210.                     c=line[i];
  211.                 if(c == '\n' || c == EOS || c == '*')
  212.                     break;            /* EOL */
  213.             
  214.                 if( (c=chartab[c]))
  215.                     seq[++(*len)]=c;
  216.                 }
  217.             if(*len == MAXLEN || c == '*') break;
  218.             }
  219.         break;
  220. /***********************************************/
  221.         case PEARSON:
  222.             while(*line != '>')
  223.                 fgets(line,MAXLINE+1,fin);
  224.             
  225.                         for(i=1;i<=strlen(line);i++)  /* DES */
  226.                 if(line[i] != ' ') break;
  227.             strncpy(sname,line+i,MAXNAMES); /* remember entryname */
  228.             sname[MAXNAMES]=EOS;
  229.             rtrim(sname);
  230.                         blank_to_(sname);
  231.  
  232.             *tit=EOS;
  233.             
  234.             *len=0;
  235.             while(fgets(line,MAXLINE+1,fin)) {
  236.                 for(i=0;*len < MAXLEN;i++) {
  237.                     c=line[i];
  238.                 if(c == '\n' || c == EOS || c == '>')
  239.                     break;            /* EOL */
  240.             
  241.                 if( (c=chartab[c]))
  242.                     seq[++(*len)]=c;
  243.             }
  244.             if(*len == MAXLEN || c == '>') break;
  245.             }
  246.         break;
  247. /**********************************************/
  248.         case GDE:
  249.             if (dnaflag) {
  250.                 while(*line != '#')
  251.                     fgets(line,MAXLINE+1,fin);
  252.             }
  253.             else {
  254.                 while(*line != '%')
  255.                     fgets(line,MAXLINE+1,fin);
  256.             }
  257.             
  258.             for (i=1;i<=MAXNAMES;i++) {
  259.                 if (line[i] == '(' || line[i] == '\n')
  260.                                   {
  261.                                     i--;
  262.                                     break;
  263.                                   }
  264.                 sname[i-1] = line[i];
  265.             }
  266.             sname[i]=EOS;
  267.             if (sname[i-1] == '(') sscanf(&line[i],"%d",&offset); /* dgg, added & to offset */
  268.             else offset = 0;
  269.             for(i=MAXNAMES-1;i > 0;i--) 
  270.                 if(isspace(sname[i])) {
  271.                     sname[i]=EOS;    
  272.                     break;
  273.                 }        
  274.                         blank_to_(sname);
  275.  
  276.  
  277.             *tit=EOS;
  278.             
  279.             *len=0;
  280.             for (i=0;i<offset;i++) seq[++(*len)] = '-';
  281.             while(fgets(line,MAXLINE+1,fin)) {
  282.             if(*line == '%' || *line == '#' || *line == '"') break;
  283.                 for(i=0;*len < MAXLEN;i++) {
  284.                     c=line[i];
  285.                 if(c == '\n' || c == EOS) 
  286.                     break;            /* EOL */
  287.             
  288.                 if( (c=chartab[c]))
  289.                     seq[++(*len)]=c;
  290.                 }
  291.             if(*len == MAXLEN) break;
  292.             }
  293.         break;
  294. /***********************************************/
  295.     }
  296.     
  297.     if(*len == MAXLEN)
  298.         warning("Sequence %s truncated to %d residues",
  299.                 sname,(pint)MAXLEN);
  300.                 
  301.     seq[*len+1]=EOS;
  302. }
  303.  
  304.  
  305. int readseqs(int first_seq) /*first_seq is the #no. of the first seq. to read */
  306. {
  307.     char line[FILENAMELEN+1];
  308.     static char seq1[MAXLEN+2],sname1[MAXNAMES+1],title[MAXTITLES+1];
  309.     int i,j,no_seqs;
  310.     static int l1;
  311.     static Boolean dnaflag1;
  312.     
  313.     if(usemenu)
  314.         getstr("Enter the name of the sequence file",line);
  315.     else
  316.         strcpy(line,seqname);
  317.     if(*line == EOS) return -1;
  318.     
  319.     if((fin=fopen(line,"r"))==NULL) {
  320.         error("Could not open sequence file %s",line);
  321.         return -1;      /* DES -1 => file not found */
  322.     }
  323.     strcpy(seqname,line);
  324.     no_seqs=0;
  325.     check_infile(&no_seqs);
  326.     printf("\nSequence format is %s\n",formatNames[seqFormat]);
  327.  
  328. /* DES DEBUG 
  329.     fprintf(stdout,"\n\n File name = %s\n\n",seqname);
  330. */
  331.     if(no_seqs == 0)
  332.         return 0;       /* return the number of seqs. (zero here)*/
  333.  
  334.     if((no_seqs + first_seq -1) > MAXN) {
  335.         error("Too many sequences. Maximum is %d",(pint)MAXN);
  336.         return 0;       /* also return zero if too many */
  337.     }
  338.  
  339. /* DES */
  340. /*    if(seqFormat == CLUSTAL) {
  341.         fprintf(stdout," \n no of sequences = %d",(pint)no_seqs);
  342.         return no_seqs;
  343.     }
  344. */
  345.  
  346.     for(i=first_seq;i<=first_seq+no_seqs-1;i++) {    /* get the seqs now*/
  347.         output_index[i] = i;    /* default output order */
  348.         if(seqFormat == CLUSTAL) 
  349.             get_clustal_seq(sname1,seq1,&l1,title,i-first_seq+1);
  350.         if(seqFormat == MSF)
  351.                 get_msf_seq(sname1,seq1,&l1,title,i-first_seq+1);
  352.         else
  353.             get_seq(sname1,seq1,&l1,title);
  354.         if(l1 > MAXLEN) {
  355.             error("Sequence too long. Maximum is %d",(pint)MAXLEN);
  356.             return 0;       /* also return zero if too many */
  357.         }
  358.         seqlen_array[i]=l1;                   /* store the length */
  359.         strcpy(names[i],sname1);              /*    "   "  name   */
  360.         strcpy(titles[i],title);              /*    "   "  title  */
  361.  
  362.         if(!explicit_dnaflag) {
  363.             dnaflag1 = check_dnaflag(seq1,l1); /* check DNA/Prot */
  364.                 if(i == 1) dnaflag = dnaflag1;
  365.         }            /* type decided by first seq*/
  366.         else
  367.             dnaflag1 = dnaflag;
  368.  
  369.         if( (!explicit_dnaflag) && (dnaflag1 != dnaflag) )
  370.             warning(
  371.     "Sequence %d [%s] appears to be of different type to sequence 1",
  372.             (pint)i,sname1);
  373.  
  374.         if(dnaflag)
  375.             n_encode(seq1,seq_array[i],l1); /* encode the sequence*/
  376.         else                    /* as ints  */
  377.             p_encode(seq1,seq_array[i],l1);
  378.     }
  379.  
  380.     fclose(fin);
  381. /*
  382.    JULIE
  383.    check sequence names are all different - otherwise phylip tree is 
  384.    confused.
  385. */
  386.     for(i=first_seq;i<=first_seq+no_seqs-1;i++) {
  387.         for(j=i+1;j<=first_seq+no_seqs-1;j++) {
  388.             if (strncmp(names[i],names[j],MAXNAMES) == 0) {
  389.                 error("Multiple sequences found with same name (first %d chars are significant)", MAXNAMES);
  390.                 return -1;
  391.             }
  392.         }
  393.     }
  394.             
  395.     return no_seqs;    /* return the number of seqs. read in this call */
  396. }
  397.  
  398.  
  399. static Boolean check_dnaflag(char *seq, int slen)
  400. /* check if DNA or Protein
  401.    The decision is based on counting all A,C,G,T,U or N. 
  402.    If >= 85% of all characters (except -) are as above => DNA  */
  403. {
  404.     int i, c, nresidues, nbases;
  405.     float ratio;
  406.     
  407.     nresidues = nbases = 0;    
  408.     for(i=1; i <= slen; i++) {
  409.         if(seq[i] != '-') {
  410.             nresidues++;
  411.             if(seq[i] == 'N')
  412.                 nbases++;
  413.             else {
  414.                 c = res_index(nucleic_acid_order, seq[i]);
  415.                 if(c >= 0)
  416.                     nbases++;
  417.             }
  418.         }
  419.     }
  420.     if( (nbases == 0) || (nresidues == 0) ) return FALSE;
  421.     ratio = (float)nbases/(float)nresidues;
  422. /* DES     fprintf(stdout,"\n nbases = %d, nresidues = %d, ratio = %f\n",
  423.         (pint)nbases,(pint)nresidues,(pint)ratio); */
  424.     if(ratio >= 0.85) 
  425.         return TRUE;
  426.     else
  427.         return FALSE;
  428. }
  429.  
  430.  
  431.  
  432. static void check_infile(int *nseqs)
  433. {
  434.     char line[MAXLINE+1];
  435.     int i;    
  436.  
  437.     *nseqs=0;
  438.     while (fgets(line,MAXLINE+1,fin) != NULL) {
  439.         if ((*line != '\n') && (*line != ' ') && (*line != '\t'))
  440.             break;
  441.     }
  442.         
  443.     for(i=0;i<=6;i++) line[i] = toupper(line[i]);
  444.  
  445.     if( linetype(line,"ID") ) {                    /* EMBL/Swiss-Prot format ? */
  446.         seqFormat=EMBLSWISS;
  447.         (*nseqs)++;
  448.     }
  449.         else if( linetype(line,"CLUSTAL") ) {
  450.         seqFormat=CLUSTAL;
  451.     }
  452.         else if( linetype(line,"PILEUP") ) {
  453.         seqFormat = MSF;
  454.     }
  455.     else if(*line == '>') {                        /* no */
  456.         seqFormat=(line[3] == ';')?PIR:PEARSON; /* distinguish PIR and Pearson */
  457.         (*nseqs)++;
  458.     }
  459.     else if((*line == '"') || (*line == '%') || (*line == '#')) {
  460.         seqFormat=GDE; /* GDE format */
  461.         if (*line == '%') {
  462.                         (*nseqs)++;
  463.             dnaflag = FALSE;
  464.             explicit_dnaflag = TRUE;
  465.         }
  466.         else if (*line == '#') {
  467.             (*nseqs)++;
  468.             dnaflag = TRUE;
  469.             explicit_dnaflag = TRUE;
  470.         }
  471.     }
  472.     else {
  473.         seqFormat=UNKNOWN;
  474.         return;
  475.     }
  476.  
  477.     while(fgets(line,MAXLINE+1,fin) != NULL) {
  478.         switch(seqFormat) {
  479.             case EMBLSWISS:
  480.                 if( linetype(line,"ID") )
  481.                     (*nseqs)++;
  482.                 break;
  483.             case PIR:
  484.             case PEARSON:
  485.                 if( *line == '>' )
  486.                     (*nseqs)++;
  487.                 break;
  488.             case GDE:
  489.                 if(( *line == '%' ) && ( dnaflag == FALSE))
  490.                     (*nseqs)++;
  491.                 else if (( *line == '#') && ( dnaflag == TRUE))
  492.                     (*nseqs)++;
  493.                 break;
  494.             case CLUSTAL:
  495.                 *nseqs = count_clustal_seqs();
  496. /* DES */             /* fprintf(stdout,"\nnseqs = %d\n",(pint)*nseqs); */
  497.                 fseek(fin,0,0);
  498.                 return;
  499.                 break;
  500.             case MSF:
  501.                 *nseqs = count_msf_seqs();
  502.                 fseek(fin,0,0);
  503.                 return;
  504.                 break;
  505.             case USER:
  506.             default:
  507.                 break;
  508.         }
  509.     }
  510.     fseek(fin,0,0);
  511. }
  512.  
  513.  
  514. static int count_clustal_seqs(void)
  515. /* count the number of sequences in a clustal alignment file */
  516. {
  517.     char line[MAXLINE+1];
  518.     int  nseqs;
  519.  
  520.     while (fgets(line,MAXLINE+1,fin) != NULL) {
  521.         if(!blankline(line)) break;        /* Look for next non- */
  522.     }                        /* blank line */
  523.     nseqs = 1;
  524.  
  525.     while (fgets(line,MAXLINE+1,fin) != NULL) {
  526.         if(blankline(line)) return nseqs;
  527.         nseqs++;
  528.     }
  529.  
  530.     return 0;    /* if you got to here-funny format/no seqs.*/
  531. }
  532.  
  533. static int count_msf_seqs(void)
  534. {
  535. /* count the number of sequences in a PILEUP alignment file */
  536.  
  537.     char line[MAXLINE+1];
  538.     int  nseqs;
  539.  
  540.     while (fgets(line,MAXLINE+1,fin) != NULL) {
  541.         if(linetype(line,"//")) break;
  542.     }
  543.  
  544.     while (fgets(line,MAXLINE+1,fin) != NULL) {
  545.         if(!blankline(line)) break;        /* Look for next non- */
  546.     }                        /* blank line */
  547.     nseqs = 1;
  548.  
  549.     while (fgets(line,MAXLINE+1,fin) != NULL) {
  550.         if(blankline(line)) return nseqs;
  551.         nseqs++;
  552.     }
  553.  
  554.     return 0;    /* if you got to here-funny format/no seqs.*/
  555. }
  556.  
  557. static void p_encode(char *seq, char *naseq, int l)
  558. {                /* code seq as ints .. use gap_pos2 for gap */
  559.     register int i;
  560. /*    static char *aacids="CSTPAGNDEQHRKMILVFYW";*/
  561.     
  562.     for(i=1;i<=l;i++)
  563.         if(seq[i] == '-')
  564.             naseq[i] = gap_pos2;
  565.         else
  566.             naseq[i] = res_index(amino_acid_codes,seq[i]);
  567.     naseq[i] = -3;
  568. }
  569.  
  570. static void n_encode(char *seq,char *naseq,int l)
  571. {                /* code seq as ints .. use gap_pos2 for gap */
  572.     register int i,c;
  573. /*    static char *nucs="ACGTU";    */
  574.     
  575.     for(i=1;i<=l;i++) {
  576.         if(seq[i] == '-')                 /* if a gap character -> code = gap_pos2 */
  577.             naseq[i] = gap_pos2;   /* this is the code for a gap in */
  578.         else {                     /* the input files */
  579.             c=res_index(nucleic_acid_order,seq[i]);
  580.             if (c == 4) c=3;    /* T = U */
  581.             if(c >= 0) naseq[i]=c;
  582.             else       naseq[i]=5;    /* N */
  583.         }
  584.     }
  585.     naseq[i] = -3;
  586. }
  587.  
  588. static int res_index(char *t,char c)
  589. {
  590.     register int i;
  591.     
  592.     for(i=0;t[i] && t[i] != c;i++)
  593.         ;
  594.     if(t[i]) return(i);
  595.     else return -1;
  596. }
  597.